Instability and Spatiotemporal Dynamics of Alternans in Paced Cardiac Tissue 
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We derive an equation that governs the spatiotemporal dynamics of small amplitude alternans in 
paced cardiac tissue. We show that a pattern-forming linear instability leads to the spontaneous 
formation of stationary or traveling waves whose nodes divide the tissue into regions with opposite 
phase of oscillation of action potential duration. This instability is important because it creates 
dynamically an heterogeneous electrical substrate for inducing fibrillation if the tissue size exceeds a 
fraction of the pattern wavelength. We compute this wavelength analytically as a function of three 
basic length scales characterizing dispersion and inter-cellular electrical coupling. 



PACS numbers: 87.19.Hh, 05.45.-a, 05.45. Gg, 89.75.-k 



It is well established that the duration of cardiac ex- 
citation can oscillate from beat to beat at sufficiently 
short pacing interval Pioneering studies by Nolasco 
and Dahlen Q and Guevara et al. || have demonstrated 
that the generic sequence LSLS... of long and short ac- 
tion potential duration (APD), known as alternans, is a 
direct consequence of the restitution relationship 



APD n+i = 



(1) 



between the APD generated by the n th + 1 stimulus, 
APD n+1 , and the diastolic time interval DP 1 during 
which the tissue recovers its resting properties after the 
end of the previous (n th ) action potential. If we denote 
the interval between the n th and n th + 1 stimulus by 
T™, we must have DP 1 = T n - APD n . Then, for a 
fixed period: T n = r for all n, Eq. [I] yields the map 
APD n = /(t - APD' 1 - 1 ) that has a period doubling in- 
stability if the slope /' of the restitution curve evaluated 
at its fixed point exceeds unity, which generically occurs 
as the period is decreased. 

Over the last decade, the study of alternans |pj— |Tl|j , and 
their control p2| , has become a main focus of research 
because of the potentially crucial link of this dynamical 
instability with cardiac fibrillation p3| . However, there 
is presently no simple analytical understanding of how 
the bifurcation to alternans is manifested spatiotempo- 
rally in paced cardiac tissue. Analytical progress to date 
is limited to the one-dimensional circulation of electrical 
impulse in a ring of tissue . 

In this letter, we derive an equation that governs the 
spatiotemporal dynamics of alternans close to the onset 
of instability. This enables us to obtain a quantitative 
analytical understanding of the formation of recently ob- 
served complex patterns of APD oscillations that can 
promote fibrillation p~pT]. A crucial feature of these 
patterns is that the APD oscillates with opposite phases 
in two (or more) spatially extended regions of tissue, i.e. 
with a sequence LSLS... in one region and SLSL... in the 
other. These "discordant alternans" have been observed 
experimentally in both two-dimensional [g and linear 
strands |J of cardiac tissue, as well as in ionic model 
simulations |p|-pl|. Moreover, they have been shown to 



lead to the formation of conduction blocks |9) as well as to 
the onset of spiral wave formation and fibrillation [|| . We 
show here that discordant alternans result from a pattern- 
forming linear instability that has interesting similari- 
ties with classic instabilities leading to the spontaneous 
formation of spatially periodic patterns in nature (such 
as Rayleigh-Benard convection, Taylor-Couette flow, etc 
p4j), but also presents some unique features. 

We consider a one-dimensional (1-d) homogeneous ca- 
ble of length L paced at period r from one end (x = 0). 
Close to the onset of instability, we can expand the APD 
and the period in the form 



APD n (x) 
T n (x) 



APD C + a(x, t) e i7r 
t c - 5t + b(x, t) e l 



(2) 
(3) 



where APD C and r c are the APD and the period evalu- 
ated at the bifurcation point of the map (/' = 1), x mea- 
sures the position along the cable, and St = t c — r <C t c . 
In this range of period, a and b vary slowly from beat 
to beat, which allows us to treat the time, t = nr, as a 
continuous variable; the fast beat-to-beat oscillations are 
contained in the exponential factor e l7rn . 

A relation between a and b can first be derived by 
noting that T n (x) is the difference of arrival time of two 
subsequent action potentials at x flllfl , or 



T n (x) =t+ r 

Jo 



dx' 



c(DI n (x')) J cpJ™- 1 ^'))' 



(4) 



where the first (second) integral on the right-hand-side 
is the time required for the leading front of the action 
potential to travel from the paced end to x at the n th 
(n th — 1) stimulus; concomitantly, c(DI) is the standard 
dispersion curve that relates the propagation speed of 
this front with the local diastolic interval. The disper- 
sion curve is typically steeply increasing at small DI and 
flat at large DI. Substituting Eqs. |^-|| in Eq. || with 
DI n (x) = T n (x) — APD n (x), and expanding to linear 
order in a and b, we obtain b(x) ~ J a(x')dx' / A where 
we have defined A = c 2 /(2c'), with c and c' = dc/dDI 
evaluated at the bifurcation, and assumed that A is much 
larger than the scale over which a varies. 
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Next, in order to derive an evolution equation for the 
amplitude a(x, t), we first neglect the influence of the 
electrical coupling between cells on the APD. This allows 
us to assume that the restitution relationship ([[]), and 
hence the second iteration of the map 



A p D n+2 = _ f(jn _ APL> n^ 



(5) 



continues to be valid even when the APD is non-spatially 
uniform; we shall soon see why it is crucial to relax this 
assumption. We substitute Eqs. in Eq. ||, and ex- 
pand the right-hand-side keeping only the dominant lin- 
ear and weakly nonlinear terms. Furthermore, we use the 
aforementioned fact that, close to onset, a varies slowly 
from beat to beat, and we therefore expand the left- 
hand-side as APD n+2 ~ APD 71 + 2da/dne l7rn , where 
da/dn = rda/dt. Finally, we use the integral relation 
between a and b derived earlier. After equating both 
sides of Eq. ||[ we obtain 



rd t a 



era 



5 a 



A 



a(x') 



(6) 



where a = f"(r - r c )/2, g = /" 2 /4 - /'"/6, and all 
derivatives are evaluated at the bifurcation point. 

In order to test this evolution equation, we simulate 
the standard cable equation 



a t v = D%v-(i bm + i ext )/c„ 



(7) 



for the membrane current I lon given by the Noble model 
Jl5| with time in units of millisecond (ms), D — 2.5 x 10~ 4 
enr/ms, C m — 12^F/cm 2 , dx=0.01 cm, dt=0.05 ms, 
and I ext modeling a sequence of stimuli applied at x = 
at the pacing interval r. To determine the parameters of 
the amplitude equation, we compute the restitution and 
dispersion curves by pacing Eq. [t] in a short cable and by 
using two subsequent stimuli spaced by different intervals 
to vary DI; we also use V = —40 mV as threshold of the 
transmembrane voltage to define the APD. We impose 
zero gradient boundary conditions on V and a at the two 
ends of the cable in all our simulations. 

Fig. |l](a) shows that Eq. ^ produces discordant al- 
ternans consistent with the picture that restitution and 
dispersion suffice to produce this state |p|-pT[. However, 
the magnitude of the spatial gradient of a at the node in- 
creases with time, and can be shown to diverge in a finite 
time. This yields an unphysical spatial discontinuity of 
APD, which also occurs if the system of coupled maps 
(Eqs. and |), fr om which the amplitude equation is 
derived, is solved numerically. Note that the analogous 
system of coupled maps for a I-d pulse circulating in a 
ring produces a smooth APD modulation starting from 
a smooth initial condition || , which highlights the much 
stronger role of dispersion at producing heterogeneity of 
APD during pacing than circulation. 

This discontinuity, which is absent in the simulation of 
the cable equation in Fig. |l](b), can be cured by adding 
spatial derivative terms to the amplitude equation. Since 
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FIG. 1. Amplitude a of APD oscillation vs x for Noble 
parameters and r = 258 ms. (a): profiles obtained with Eq. 
^| at different times (dotted lines) and final stationary profile 
(solid line), (b): stationary profiles obtained with Eq. (7 (solid 
line) and Eq. H (dashed line). Nodes (a — 0) separate tissue 
regions with n out of phase oscillations. 



the cable is paced at one end, the underlying basic state 
(i.e. traveling pulses) is not invariant under parity sym- 
metry. Hence, in addition to d 2 a, a term proportional 
to d x a must generally be included, which yields the final 
form of our amplitude equation 



rdta 



aa — ga 



x dx^_ 
A 



a(x')-wd x a + i 2 dla. (8) 



To calculate the new lengthscales w and £, we must deter- 
mine how the electrical coupling between cells modifies 
the restitution relationship (Eq. ||) . For complex electro- 
physiological models such as Noble, this generally needs 
to be done numerically using a procedure that will be dis- 
cussed elsewhere. For a simple two-variable ionic model, 
which we study below, the analytical expressions 



w = 2D/c, 
£ = (Dx APD.) 1 / 2 , 



(9) 
(10) 



can be derived by interpreting Eq. as a diffusion equa- 
tion with a source Ii on , and expressing V as space-time 
integral of GIi on , where G is the standard Green's func- 
tion of the diffusion equation. This integral can then be 
evaluated analytically because the action potential shape 
is simply triangular for this model, and used to calculate 
w and £. Eq. [w] has the simple physical interpreta- 
tion that V diffuses a length ~ £ in the time interval of 
one APD. Therefore, the repolarization of a given cell 
is influenced by other cells within a length ~ £ of ca- 
ble. In addition, repolarization of this cell is influenced 
unequally by its left and right neighboring cells because 
these cells are activated at different times by the prop- 
agating wavefront. Clearly, this asymmetry must vanish 
in the limit c — > 00 where all cells are activated simul- 
taneously consistent with Eq. ||. Fig. |l|(b) shows that 
our regularized amplitude equation (||) now produces a 
smoothly varying stationary profile of a that agrees well 
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TABLE I. Values of various lengths in cm with Xtheor (Eq. 
|l2] or Eq. |l3"| and A s ; m (from simulations of Eq. Q). 



Model 


A 


w 


t 


Ath e0 7'/4 


^sim / 4 




Noble 


49.1 


0.045 


0.18 


2.33 


2.6 


2.75 


Two- variable 


3.55 


0.031 


0.235 


1.33 


1.1 


1.15 



with the simulation of the cable-Noble equation, where a 
is obtained from the APD using Eq. |^. 

The genesis of discordant alternans can be understood 
by computing the linear stability spectrum of the spa- 
tially homogeneous state (a = 0). We have calculated 
this spectrum both numerically for different L, and ana- 
lytically for the large L limit. The main result is that the 
wave pattern can emerge from the amplification of either 
a unique finite wavelength mode, which yields a station- 
ary pattern, or from a discrete set of complex modes that 
approach a continuum in the limit L — > oo, and yields a 
traveling pattern. There is indeed experimental evidence 
for both stationary || and traveling || waves. 

We treat here explicitly the large L limit since it pro- 
vides the basis to understand fmite-L patterns. In this 
limit, we can analyze stability by differentiating Eq. |§| 
with respect to x and letting a(x, t) ~ ^kx+Qt^ w ^ n k th 



f2 = fl r + iQi and k = k r + iki 
once the eigenvalue equation 



complex, which yields at 



£ A kr -i[wk-l/(Ak)} 



(11) 



When dispersion is weak, a unique real mode with fcj = 
grows faster than the other complex modes and yields a 
stationary pattern. Its wavelength A = 2n/k r , is deter- 
mined by the condition f2j — 0, which yields 



A = 2tt(«;A) 1/2 , 



(12) 



in good agreement with the wavelength observed in sim- 
ulations of the cable-Noble equation (Table |). Note that 
cos k r x is an exact eigenvector of Eq. || linearized around 
a = that satisfies d x a = at the two cable ends when 
L is an integer multiple of A/2. The threshold of insta- 
bility occurs when fl r = 0, or for a period r t h defined by 

<Tth = f"(Tth-T c )/2 = e/(wA). 

In the opposite limit where dispersion is strong, com- 
plex modes that grow exponentially at large x are the 
most unstable. It is simple to deduce from Eq. [Il] that 
each &;-mode travels towards the pacing end of the cable, 
but a wave-packet constructed from a linear superposi- 
tion of these modes has a group velocity that makes the 
packet move away from the pacing end. This is the signa- 
ture of a convective instability [14] where perturbations 
are transported as they grow, similarly to Taylor-Couette 
vortices developing in an axial flow |16|. In such a sit- 
uation, patterns are only transient unless the group ve- 
locity vanishes, or dfli/dk r — 0, and hence they grow 
at a fixed point in space. Moreover, the fastest grow- 
ing wavelength that dominates at large time must corre- 
spond to a maximum of f2 r , which yields the additional 
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FIG. 2. Space-time plots of a obtained by simulations of 
Eq. ^| for parameters of the two- variable model, showing ab- 
solutely unstable (left and r = 295 ms) and convectively un- 
stable (right and r = 298 ms ) wave patterns. The crosses 
denote the positions of the nodes (a — 0) . 

condition dfl r /dk r = 0. These two conditions are equiv- 
alent to requiring that dQ/dk — 0. From this, we deduce 
that the threshold of absolute instability occurs when 
o~th — (3/2)(£/2A) 2 / 3 , with a pattern of wavelength 



A = (4tt/V3)(2C 2 A) 1 / 3 , 



(13) 



which travels with phase velocity f2jA/(27r) where f2j = 
(3V3/2)(£/2A) 2 / 3 . It can also be deduced that traveling 
waves are favored over a stationary pattern when A = 
2c 2 /c' <C £ 4 /ui 3 , and hence for large enough dispersion 
(c'), and vice versa in the opposite limit. 

Nonlinearities are mainly responsible for saturating the 
amplitude of the linear modes, but both the wavelength 
and evolution of the pattern can generally vary with dis- 
tance from onset. There is a close analogy between the 
amplitude equation (JsJ) and the real Ginzburg-Landau 
equation that has been extensively studied in the con- 
text of phase transitions and front propagation |l4|] . The 
dynamics is richer here because the integral term origi- 
nating from dispersion causes a non-local interaction of 
the fronts separating two out-of-phase oscillating regions 
with the pacing end of the cable. 

As a final test of our theory, we study a two-variable 
model where all the parameters of the amplitude equation 
can be calculated analytically including w and £ given by 
Eqs. | 10, and for which our theory predicts traveling 
waves. In this model, Iion/C m is the sum of a slow out- 
ward current, I„/C m = t^ 1 (S + (1 — S) V/V c ), and a fast 



inward current, Ii/C„ 
the latter is controlled by 



t„ hS, where inactivation of 



d t h = (1 - S - h)/(r-(l -S) + St+) 



(14) 



In addition, V is dimensionless and S = (1 + tanh((y — 
V c )/e))/2. We choose V c = 0.1, r = 150 ms, r a = 6 ms, 
t~ = 60 ms, t + = 12 ms, e — 0.005, and we simulate Eq. 
|?] with the same D as before, dx=0.01 cm, dt=0.02 ms, 
and V = 0.1 as threshold for the APD. 



3 



♦♦♦♦♦♦( 

♦ « 

♦ 4 

♦ < 

♦♦♦♦♦♦ 



♦♦♦♦♦♦ 



BOljOO . 

lOOOO 
)0000 
OOOOO 
000 

ooo 
ooo 
ooo 
c o o 
ooo 
ooo 
poo 
" oo 
'ooo 
** S o o o o 
po 






100 









•100 




50 


(m 







-50 




50 




-50 



x=290ms 



0.5 









x=290ms ; 




2S0 



290 300 310 

x (ms) 



32(1 



12 3 4 

x (cm) 



FIG. 3. Stability diagram of two-variable cable model with 
domains of no-alternans (open circles), concordant alternans 
(filled squares), and discordant alternans (filled diamonds); 
conduction blocks form at smaller r not shown here. Bound- 
aries between the same domains obtained by simulations of 
the amplitude equation (^) are shown by solid lines. The 
dashed line denotes the bifurcation period for alternans pre- 
dicted by the map of Eq. [[]. On the right panel we show 
profiles of a vs x at times t (solid), t + 2r (dashed) and t + 4r 
(dotted lines), for three different cable lengths. 



As predicted, we observe traveling waves in this model 
with a wavelength that agrees well with Eq. 



13. Fig. 



illustrates sustained and transient wave patterns above 
(left) and below (right) the onset of absolute stability, 
respectively. Furthermore, Fig. shows that the results 
of simulations of the cable and amplitude equations are 
in good agreement over a wide range of L and r, in- 
cluding the fact that the onset of instability occurs at a 
shorter period in a larger cable (and, hence, for a slope 
of the restitution curve larger than unity). Fig. [| is 
qualitatively similar for stationary waves, but this over- 
stabilization is smaller because dispersion is weaker. 

The APD oscillations at a fixed x are periodic when 
the wave is stationary and quasiperiodic when it travels. 
For both cases, we find here that A is independent of 
L, in contrast to the oscillations produced by a pulse 
circulating in a ring where it is known that A ~ 2L/(2i + 
1) for weak dispersion |5| (with i integer and L = ring 
perimeter). As will be discussed elsewhere, our theory 
applied to the ring shows that the bifurcation to alternans 
is finite dimensional with i = being the most unstable 
mode, in agreement with the fact that it is this mode 
that is generically selected in experiments [Q and ionic 
model simulations M nn. In addition, it shows that the 
gradient term (— wd x a) can lead to quasiperiodicity even 
in the absence of dispersion (c' = 0). 

Our results demonstrate that the formation of discor- 
dant alternans is crucially affected by the effect of electri- 
cal coupling (diffusion) on repolarization, in addition to 
restitution and dispersion. Dispersion is responsible for 
the formation of nodes and spatial gradients of APD that 
steepen with time. Diffusion, in turn, tends to spread the 
APD spatially, and also induces a drift of the pattern 
away from the pacing site that is induced by the more 
subtle gradient term (—wd x a) in the amplitude equation. 



When dispersion is sufficiently weak, drift balances dis- 
persion and produce a stationary pattern. In the opposite 
limit, the tendency for dispersion to form steep gradients 
of APD is balanced by the spreading effect of diffusion. 
Nodes then travel, cyclically disappearing (appearing) at 
the pacing (opposite) end of the cable. 

In conclusion, we have derived a simple evolution equa- 
tion that describes the universal spatiotemporal dynam- 
ics of small amplitude alternans in paced cardiac tissue. 
Moreover, we have shown that discordant wave patterns 
|p[-pT| , which are linked to fibrillation || , result from a fi- 
nite wavelength linear instability. Hence, their formation 
requires a minimum tissue size L m i n ~ A/4, required for 



at least one node to form. The value of L r 



that we 



measure in simulations of reaction-diffusion models are 
actually close to A/4 with A predicted by Eq. [f2] and 
Eq. 13, respectively (Table ||). This equation can be ex- 
tended to higher dimensions as well as to include hetero- 
geneities. Our preliminary results in a 2-d homogeneous 
tissue are very similar to 1-d. Finally, this equation is 
readily applicable to model a non-constant pacing inter- 
val and provides a theoretical basis to study the control 
of alternans in spatially extended tissue. 
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